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Abstract 

The unprecedentedly bright afterglow of Swift GRB 130606A at z = 5.91 gave us a unique 
opportunity to probe the reionization era by high precision analyses of the redward damping 
wing of Lya absorption, but the reported constraints on the neutral hydrogen fraction (/ H i) in 
intergalactic medium (IGM) derived from spectra taken by different telescopes are in contra¬ 
diction. Here we examine the origin of this discrepancy by analyzing the spectrum taken by 
VLT with our own analysis code previously used to fit the Subaru spectrum. Though the VLT 
team reported no evidence for IGM H I using the VLT spectrum, we confirmed our previous 
result of preferring non-zero IGM H I (the best-fit /hi ~ 0.06, when IGM H I extends to the GRB 
redshift). The fit residuals of the VLT spectrum by the model without IGM HI show the same 
systematic trend as the Subaru spectrum. We consider that the likely origin of the discrepancy 
between the two teams is the difference of the wavelength ranges adopted in the fittings; our 
wavelength range is wider than that of the VLT team, and also we avoided the shortest wave¬ 
length range of deep Lya absorption (A 0 bs < 8426 A), because this region is dominated by H I in 
the host galaxy and the systematic uncertainty about host H i velocity distribution is large. We 
also study the sensitivity of these results to the adopted Lya cross section formulae, ranging 
from the classical Lorentzian function to the most recent one taking into account fully quantum 
mechanical scattering. It is found that the preference for non-zero IGM Hi is robust against 
the choice of the cross section formulae, but it is quantitatively not negligible and hence one 
should be careful in future analyses. 

Key words: Techniques: spectroscopic — Gamma-ray burst: individual: GRB 130606A— dark ages, 
reionization, first stars 


1 Introduction 

Hydrogen gas in the intergalactic medium (IGM) was reionized 
at redshifts z> 6, and it is important to reveal when and how 


this process occurred quantitatively for better understanding of 
the formation and evolution of the earliest galaxy populations 
(Barkana & Loeb 2007; Meiksin 2009; Robertson et al. 2010; 
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Fan 2012 for reviews). Gamma-ray bursts (GRBs) are an im¬ 
portant probe for the distant universe reaching the epoch of cos¬ 
mic reionization (z > 6), giving us useful information such as 
the cosmic star formation rate (see Piran 2005; Meszaros 2006; 
Zhang 2007; Gehrels et al. 2009 for reviews). When a high 
precision optical/near-infrared spectrum is taken for a GRB af¬ 
terglow in the reionization era, light at wavelengths shorter than 
the redshifted Ly a is almost completely absorbed by neutral 
hydrogen in IGM, and so-called Gunn-Peterson troughs ap¬ 
pear. Because of the strength of Lya absorption, GP troughs 
give only a weak lower bound for the IGM neutral fraction as 
/hi = nm/nu > 10 -3 . However, the shape of redward Lya 
damping wing can be used to further constrain /hi and hence 
the reionization history, because IGM HI would affect the shape 
of the wing, in addition to the wing by HI in the host galaxy 
(Miralda-Escude 1998). GRBs have a few advantages against 
quasars about this test. They are expected to occur in more nor¬ 
mal and less dense regions than quasars, and short duration of 
GRB emission does not affect the ionization status around host 
galaxies. The simple power-law spectrum of GRB afterglows is 
suitable for a precise fit to the damping wing shape. 

However, use of GRBs for reionization study has been ham¬ 
pered by the low event rate of high-s GRBs that are bright 
enough for high precision damping wing analyses. A weak 
upper bound on IGM Hi was obtained for GRB 050904 at 
2 = 6.3 (Kawai et al. 2006; Totani et al. 2006), but meaning¬ 
ful constraints on reionization could not be derived by GRBs at 
even higher redshifts (Greiner et al. 2009; Patel et al. 2010; 
Salvaterra et al. 2009; Tanvir et al. 2009; Cucchiara et al. 
2011) because of their low signal-to-noise ratio. The discovery 
of GRB 130606A (Ukwatta et al. 2013; Golenetskii et al. 2013; 
Castro-Tirado et al. 2013) at z = 5.91 provided us with the best 
opportunity so far for this purpose, by its exceptional bright¬ 
ness and relatively low HI column density in the host galaxy 1 . 
However, the reported constraints on /hi using spectra taken by 
different telescopes are controversial. Chornock et al. (2013) 
derived an upper bound of /hi < 0.11 (2<r) using the spectrum 
taken by Gemini. In contrast, Totani et al. (2014, hereafter 
Paper I) found a ~ 3a evidence for IGM Hi with /hi > 0.08 
from the spectrum taken by Subaru. It was argued that the 
choice of the baseline power-law index /3 = —1.99 (/„ oc 
in Chornock et al. is not supported by the observed optical/NIR 
colors favoring /3 ~ — 1. This has been confirmed by the optical- 
NIR spectrum by VLT, reporting /3 = —1.02 ± 0.03 (Hartoog et 
al. 2014, hereafter HI4). 

However, this is not the end of the story. The analysis of the 
damping wing by H14 using the VLT spectrum found no evi¬ 
dence for IGM Hi, setting an upper limit on /hi < 0.03 (3<r). 

1 After this event, GRB 140515A was detected at z = 6.33 and constraints 
on reionization have been derived (Chornock et al. 2014; Melandri et al. 
2015), but the afterglow was not as bright as that of GRB 130606A. 


This is statistically inconsistent with the result of our Paper I, in¬ 
dicating the systematic uncertainties in the Subaru/VLT spectra 
and/or analysis methods. GRB 130606A has proven that a high 
precision damping wing analysis for GRB spectra in the reion¬ 
ization era is indeed possible, and understanding systematics in 
such analyses is crucial for future studies using more GRBs at 
higher redshifts to derive reliable constraints on reionization. 
Therefore we decided to analyze the VLT spectrum by our own 
code used for the fitting to the Subaru spectrum in Paper I, to 
reveal the origin of the discrepant results. This is the primary 
aim of this work. 

Another aim of this work is to examine the effect of adopted 
formulae for the Lya cross section as a function of wavelength. 
Several different formulae have been used in the literature to 
calculate the damping wing shape by HI in a host galaxy and 
IGM, including the simplest Lorentzian (or the Voigt profile 
when convolved with the Gaussian velocity distribution), and 
the two-level approximation formula by Peebles (1993, here¬ 
after P93). In high precision damping wing analyses, difference 
of the adopted formulae may result in systematic biases. Here 
we repeat our analyses using different cross section formulae, 
including the most recent formula by Bach & Lee (2015, here¬ 
after BL15) that takes into account the fully quantum mechani¬ 
cal scattering based on the second-order time-dependent pertur¬ 
bation theory, and see how the results change. 

Unless otherwise stated, we adopt the same model parameter 
values as those in Paper I. The HI column density in the host 
galaxy, N^° st , is expressed in units of cm -2 . 

2 On the Controversial Results between 
Subaru and VLT 

2.1 Comparison between the Subaru and VLT 
Spectra 

First we directly compare the Subaru and VLT spectra of GRB 
130606A, which were taken during 10.4-13.2 and 7.2-8.7 hr 
after the burst, respectively. We use the one-dimensional VLT 
spectrum and its statistical error, which were reduced by the 
VLT team and provided to us. The original VLT spectrum has 
a wavelength binning size of 0.20 A, which is smaller than 0.74 
A for the Subaru spectrum. To compare the two spectra on the 
same wavelength grids, the VLT spectrum was converted into 
that on the Subaru grids. This was done first by associating 
every VLT grid to the closest Subaru grid, and then by taking the 
average over the VLT grids associated to a Subaru grid. The two 
spectra on the same grids are then shown in Fig. 1. The la error 
at regions without strong airglow emission is typically 1.0% and 
1.5% of the continuum level for Subaru and VLT, respectively, 
on this same wavelength grids. The VLT spectrum shows finer 
structure because of the better spectral resolution (R = A/A A ~ 
2200 and 8700, respectively). The ratio of the two spectra is 
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Fig. 1 . Top: the Subaru spectrum of the afterglow of GRB 130606A reported 
in Paper I. The 1 a error is shown, but multiplied by a factor of five for presen¬ 
tation purpose. Middle: the same as the top panel but for the VLT spectrum 
of H14. Here, the VLT spectrum has been converted on the wavelength grids 
of the Subaru spectrum (see text). Note that the spectral resolution of the 
VLT spectrum is better than the Subaru even though the wavelength grids 
are the same. Bottom: the flux ratio of the VLT to the Subaru spectrum. 

also shown in this figure. Though there are fluctuations on ~ 10 
A scale, no systematic trend is found on scales longer than 100 

A. 

2.2 Results of the Damping Wing Fits 

In Paper I, we considered two different models for the IGM HI 
distribution; in one model, the IGM HI is assumed to extend up 
to the same redshift as the GRB, i.e., zigm,u = zgrb = 5.9131, 
while in the other model, zigm,u is set at 5.83 motivated by 
a deep GP trough observed at z = 5.67-5.83 to the sightline 
of this GRB (Chornock et al. 2013). 2 It should be noted here 
that the damping wing shape is calculated assuming that IGM 

2 The fit is insensitive to the redshift lower bound zigm.i if zigm.i -C 
zigm,u, and contribution to the damping wing is mostly from a redshift 


HI is uniformly distributed with a constant comoving density, 
which is obviously too simple compared with the expected re¬ 
alistic distribution. Though it is not the purpose of this pa¬ 
per to examine this effect, this should be kept in mind when 
the implications for reionization are discussed. The best-fit 
model parameters are (/hi,/ 3) = (0.086jlo;oii,— 0.931 q; and 
(0.471 q;q®, —0.74to)o?) for the former and latter, respectively. 
The x 2 difference between the two models is not statistically 
significant, but both the two models show a ~ 3cr statistical 
preference to the model without IGM HI. It should be noted 
that some other possibilities of reddening absorption [extinction 
in the host galaxy or intervening damped Lya systems (DLAs)] 
have been ruled out in Paper I. 

However, H14 reported /? = —1.02 ± 0.03 from their op- 
tical+NIR spectrum, indicating that the low zigm.u model in 
Paper I is now disfavored. Therefore in this paper we consider 
only the zigm.u = zgrb model with the power index fixed at 
/3 = —1.02, to compare the Subaru and VLT spectra 3 . Then 
the model parameters are: the column density of HI in the host 
galaxy A r H° st , the la Gaussian velocity dispersion cr„ of Hi in 
the host galaxy, /h i of IGM H I, and the overall normalization 
factor. The fit results of the host Hi only model (/hi fixed to 
zero) and the host+IGM HI model (/hi treated as a free param¬ 
eter) to the Subaru and VLT spectra are shown in Table 1 and 
Figs 2 and 3. Here, we used the VLT spectrum converted onto 
the wavelength grids of the Subaru spectrum as described in the 
previous section, and exactly the same analysis procedures as 
Paper I were adopted to the two spectra. 

The fit results to the Subaru spectrum with the fixed value of 
/3 = —1.02 are similar to the zigm.u = zgrb model with free 
/3 reported in Paper I; a non-zero IGM contribution with /hi ~ 
0.06 is favored than the host HI only model. The relative flux 
difference between the two models is only ~ 0.5% level, but the 
overall statistical significance is (A^ 2 ) 1 ^ 2 = 3.8 ct thanks to the 
small statistical error of the spectrum. It should be noted that the 
errors on the Subaru spectrum are uncorrelated among different 
wavelength bins and the scatter is consistent with the Gaussian 
(Paper I). There are 68 data points in the fit, and there are four 
model parameters in the host+IGM HI model. The \ 2 value of 
80.62 for 68 — 4 = 64 degrees of freedom is not unreasonable 
(8% chance probability of getting a larger value). 

The reason why the IGM HI component improves the fit can 
be seen in the residuals (normalized by the observed flux error) 
shown in Fig. 2. The fit residuals of the host only model tend 
to be positive at longer wavelength range of 8800-8900 A, but 
tend to be negative at 8480-8560 A, indicating that another ab¬ 
sorption component is necessary to make the spectrum further 

interval of Az ~ 0.1 from zi G m,u [see Fig. 7 of Totani et al. (2006)], 
corresponding to a proper distance of ~ 6 Mpc. 

3 The error of 0 reported by HI 4 is smaller than the typical statistical uncer¬ 
tainties of free 0 fits in Paper I, and hence we do not marginalize the error 
of the HI 4 0 measurement. 
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Table 1 . The best fit parameters of the fittings to the Subaru and VLT spectra* 


model 

logic 

< 7 v (km/s)-I- 

IGM / H , 

2§ 

x 2 

(N 

* 

<1 

fit to the Subaru spectrum 

host HI only 

1 Q 077+O.OO8 
/ _ 0.015 

0 o+ 89 - 9 
u -0.0 

fixed to zero 

95.10 

14.48 

host+IGM Hi 

19 768~*~°' 032 
/oo _ 0 032 

62 n - *" 38 0 

OZ -U_ 62.0 

0 061” 1 ” 0 ’" 7 
u.uoi_o 007 

80.62 

- 

fit to the VLT spectrum 

host HI only 

19 806 +U U14 

i?.ouu_ 0 016 

o.o^ 0 u 

fixed to zero 

292.57 

11.89 

host+IGM Hi 

19 621“*" 0 059 

i?.ozi_ 0 Q57 

0.0 + IT 

U.Uo /Q 029 

280.68 

- 


*The fit results for the Lya damping wing models including only HI in the host galaxy (“host H1 only”) 
and including host plus IGM H 1 (“host+IGM Hi", with IGM HI extending up to the redshift of the 
GRB). The quoted errors are statistical \a. 

tThe neutral hydrogen column density in the GRB host (fV^° st ) in units of cm -2 . 

^The survey range of cr v (lrr of the Gaussian velocity distribution of host Hi) is limited to 0-100 km/s, 
motivated from the observed widths of the metal absorption lines (Paper I). 

''' The number of data points is 68. 

^ The \ 2 excess of the host H1 only model against the host+IGM HI model. 


redder. It is impossible to do this by increasing the amount of 
host H I, because it would result in too large absorption in the 
shortest wavelength range of 8420-8460 A. Since HI in IGM 
is more distant from the GRB than HI in the host, the damping 
wing by IGM HI has weaker wavelength dependence than that 
by host HI (see Fig. 2 of Paper I). This is the key for the fit im¬ 
provement by IGM H I. It should also be noted that a reduction 
of the scatter of residuals by IGM HI is seen in 8420-8460 A. 
This indicates that the damping wing shape with IGM HI is in 
better agreement with the observed data than the host HI only 
model, giving a further support for the model with IGM H I. 

Now we examine the fit results to the VLT spectrum. We 
confirmed the same result as that to the Subaru spectrum; non¬ 
zero IGM HI contribution is favored with the best fit value of 
/hi = 0.0871Q Q29 which is statistically consistent with /hi = 
0.0611q found with the Subaru spectrum. The residuals of 
the fit to the VLT spectrum in Fig. 3 show the same system¬ 
atic trend as those of the Subaru fit (in Fig. 2). Therefore 
the evidence for IGM HI contribution to the observed damp¬ 
ing wing reported in Paper I is further strengthened by an in¬ 
dependent spectrum taken by VLT. The x 2 difference implies 
that the host+IGM HI model is preferred to the host only model 
with a statistical significance of 11.89 1 / 2 = 3.4<r. However, the 
residuals shown in Fig. 2 are on average significantly larger that 
the expectation of the Gaussian distribution, suggesting that the 
errors calculated by the VLT analysis pipeline are underesti¬ 
mated. We discussed about this with the VLT team, and it is 
most likely because the lower signal-to-noise ratio with respect 
to the Subaru spectrum makes the systematic noise due to the 
bright sky emission lines more prominent. The VLT spectrum 
was taken with an exposure time considerably shorter than that 
for Subaru, and the air mass was also larger. The median of 
absolute values of the residuals for the host+IGM HI model is 
1.221, while 0.674 is expected for the Gaussian distribution. If 
we scale the error size by the ratio of 1.221/0.674, the signif¬ 
icance is reduced to 1.9cr. In this work we do not discuss the 


statistical significance about the VLT spectrum further. 

These results indicate that the discrepant conclusions (evi¬ 
dence for IGM HI by Paper I, while no evidence found by HI4) 
are a result of difference in analysis methods, rather than the 
systematic difference between the two spectra. Concerning this, 
we note that the wavelength range fitted by H14, 8406-8462 A 
(corresponding to relative velocities of 0-2000 km/s to Lya), 
is significantly shorter compared with 8426-8900 A adopted in 
our Paper I. It should be noted that the evidence for IGM HI in 
our analysis is found in the spectral shape spanning in the whole 
wavelength range. Another difference is that we excluded the 
shortest wavelength range of A 0 bs < 8426 A. As shown in Fig. 
2 of Paper I, the HI absorption in this region is dominated by 
the host HI rather than the IGM component with /h i ~ 0.1. The 
observed spectrum rapidly drops with decreasing wavelength at 
A 0 bs < 8426 A, making a precise model fitting difficult because 
of uncertainties about the velocity distribution of HI in the host 
and spectral resolution. Though the Gaussian velocity distri¬ 
bution has been taken into account for the host H I, such a pure 
Gaussian distribution is unlikely in a realistic galaxy, as inferred 
from high resolution spectroscopic studies of DLAs (Wolfe et 
al. 2005; see also Paper I for more discussion). This is why 
in Paper I we used only wavelength longer than 8426 A; this 
boundary was determined so that the fit result becomes mostly 
insensitive to the host Hi velocity dispersion, cr v (see Table 1). 

To test the effect of the wavelength range used, we per¬ 
formed the damping wing fit to the Subaru spectrum with the 
same wavelength range as H14. Indeed, we found a similar 
result to H14; the best fit is /hi = 0 with upper bounds of 
/hi < 0.0034 and 0.088 at 1 and 2cr, respectively, log 10 N^° at 
= 19.855±0.009, and cr„ = 61.8 ± 3.3 km/s. The small error of 
o v indicates that the fit in this region is highly sensitive to the 
host HI velocity distribution, but systematic uncertainty should 
be large if the distribution is not a simple Gaussian. 
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wavelength [A] 

Fig. 2. Top: the observed spectrum taken by Subaru is shown in thin solid line. The best-fit curves of the host H I only model and the host+IGM H I model are 
shown by thick dashed and solid curves, respectively. Bottom: the fit residuals (/ oba - /model)/fobs are shown for the two models. The gray shaded regions 
indicate the wavelength ranges removed from the fits because of apparent features of absorption lines or airglow. 


3 Testing Different Lya Cross Section 
Formulae 

In Paper I we used the two-level approximation formula of P93 
for the Ly a cross section as a function of wavelength. In the cal¬ 
culation of the damping wing by IGM H I, we further used the 
analytic formula of Miralda-Escude (1998), which also assumes 
the P93 formula and uses an approximation valid when the fre¬ 
quency is far from the resonance (Air\v — u a \ r a ) to make 
the analytic integration possible, where r a = 6.262 x 10 s s -1 
is the damping constant of Lya. Though this approximation is 
very good (corresponding to |A — A a |/A a 2.02 x 10 -8 ) for 
the practical damping wing analysis, we revised our code to nu¬ 
merically integrate the IGM damping wing, to calculate it for 
arbitrary formulae of Lya cross section. 

Here we test three different formulae of Lya cross section: 
the Lorentzian, the classical Rayleigh scattering formula, and 
the fitting formula of BL15 taking into account the quantum 
mechanical scattering effect, in addition to the P93 formula used 


in Paper I. (See BL15 for the exact expressions of these formu¬ 
lae.) Fig. 4 shows the change of HI optical depth as a function 
of wavelength for different formulae, using the best-fit parame¬ 
ters of the host+IGM HI model to the Subaru spectrum of GRB 
130606A. The change relative to the Lorentzian becomes larger 
with increasing wavelength, reaching 10-20% at 8900 A. The 
Rayleigh and P93 formulae result in smaller optical depth, but 
the BL15 formula larger, compared with the Lorentzian. Then 
we repeated our damping wing fitting analysis on the Subaru 
spectrum for the host HI only and host+IGM HI models, and 
X 2 are presented in Table 2. The best-fit value of /hi changes 
by less than 8%. The BL15 formula results in the smallest sig¬ 
nificance of the statistical preference for the non-zero /hi to the 
host Hi only model, but still it is larger than 3.1 ct. 

4 Discussion and Conclusions 

This paper investigated the origin of the discrepant results about 
the constraint on IGM Hi fraction from Lya damping wing 
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wavelength [A] 

Fig. 3. The same as Fig. 2, but for the VL.T spectrum converted onto the wavelength grids of the Subaru spectrum. 


Table 2. The best fit parameters in the fittings with different Ly a cross section 
formulae* _ 


cross section formula 

logio^r 4 ) 

a v (km/s) 

IGM/h, 

x 2 

Ax 2 

host HI only model 

Lorentzian 

19.8691qqJq 

0 n+ 70 - 2 

u - u -0.0 

fixed to zero 

91.81 

10.74 

Rayleigh 

19.875+°;°!° 

99 1 +03.1 
z ' z " l -22.1 

fixed to zero 

94.21 

13.50 

Peebles 

19.8771°;°“ 

0 o+ 89 - 9 

U -0.0 

fixed to zero 

95.10 

14.48 

Bach & Lee 

19-866l°;°09 

q q+63.5 

u -o.o 

fixed to zero 

90.66 

9.88 

host + IGM HI model 

Lorentzian 

19 755~*~ u ' u3a 
17 , ' jj _ 0.0 33 

100 0 +uu 
iuu.u_ 100 q 

0 057+O-9O12 

u.uj /_ 0 007 

81.07 

- 

Rayleigh 

19 765“*" 0 033 
17,/uj _0.033 

54 6+ 45 ' 4 

54.6 

0 060 + ° 008 
u.uou_ 0 007 

80.71 

- 

Peebles 

19 768~*” 0 ' 032 
/O5_o.o32 

62 o+ 38 -° 
OZ.O-62.0 

0 061 +0 ‘ 007 
U.UOl _q 007 

80.62 

- 

Bach & Lee 

19 7S1 +0-029 

100 0+0-0 
iuu.u 100 o 

0 056” 1 " 0,011 
U.UjO—0,006 

80.78 

- 


* The fitting results of the host HI only and host+IGM HI models for the Lya damping wing, using the Subaru 
spectrum. See also Table 1 for more explanations. 
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Fig. 4. Top: the flux attenuation factor e -T by the optical depth r of absorp¬ 
tion by Hi, separately for the host galaxy Hi and IGM Hi with the best-fit 
parameters of the host+IGM HI model to the Subaru spectrum (Table 1), us¬ 
ing the Lorentzian cross section formula. Middle: the fractional difference of 
t by H i in the host galaxy when various Lya cross section formulae (o-p: 
Peebles, <jr : Rayleigh, a B l - Bach-Lee) are used instead of the Lorentzian 
(ox). Bottom: the same as the middle panel, but for r by IGM Hi. 


analyses of the Subaru and VLT spectra for GRB 130606A at 
2 = 5.913. Direct comparison between the two spectra does 
not show any systematic difference on the overall shape in the 
wavelength range of 8426-8900 A used in Paper I. We repeated 
exactly the same analysis as Paper I for the VLT spectrum con¬ 
verted onto the wavelength grids points of the Subaru spec¬ 
trum, and confirmed our previous result of more than 3<r sta¬ 
tistical preference for non-zero contribution from IGM HI with 
/hi ~0.05-0.1 (assuming uniform distribution of IGM Hi up 
to the redshift of the GRB). The absorption by IGM is dom¬ 
inated by Hi within A 2 ~ 0.1 from the upper redshift bound 
(“iGM,u = zgr.b), corresponding to a proper distance of ~ 6 
Mpc. The host H 1 only model is disfavored because of the sys¬ 
tematic trend of fit residuals (positive at 8640-8900 A, while 


negative at 8450-8560 A) indicating that another absorption 
component is necessary to make the spectrum redder, in addi¬ 
tion to Hi in the host galaxy. This trend has been confirmed 
also using the VLT spectrum. 

Therefore we consider that the discrepant results were ob¬ 
tained because of different analyses methods. H14 adopted the 
wavelength range of 8406-8462 A, which overlaps only with 
the shortest part of the range adopted by us in Paper I, and H14 
also included the deep Lya absorption region down to the res¬ 
onant Ly a wavelength (8406 A). In Paper I we excluded the 
wavelengths shorter than 8426 A because this region is highly 
sensitive to the velocity distribution of H 1 in the host galaxy, 
which is rather uncertain and unlikely to be a pure Gaussian. 
Indeed, we found a result consistent with H14 (i.e., no evidence 
for non-zero /hi) when we adopted the same wavelength range 
as H14. 

Finally, we examined the robustness of these results against 
the adopted Lya cross section formulae, by testing the 
Lorentzian, Rayleigh scattering formula, and the latest BL15 
fitting formula taking into account a fully quantum mechanical 
scattering effect, in addition to the Peebles’ two-level approx¬ 
imation formula used in Paper I. We found that the preference 
for the non-zero /hi is robust against the difference of these 
formulae. The BL15 formula results in the lowest statistical 
significance, but still it is greater than 3.1 <7 and the best-fit /hi 
changes at most 8 %. However this effect is quantitatively not 
negligible, and one must be careful about the Lya cross section 
formulae in future studies. 

These results demonstrate that high precision Lya damp¬ 
ing wing analyses of high -2 GRB afterglows are a powerful 
approach to measure IGM H 1 fraction, which is sensitive to a 
small value of /hi ~ 0.05. The same results were obtained by 
the two spectra taken by the two different telescopes, when the 
model fitting method is the same, indicating that the system¬ 
atic uncertainty of the spectral measurement can be controlled 
to allow discrimination of sub-percent level relative flux change 
in the spectral shape. However, a high precision fitting is eas¬ 
ily affected by systematic uncertainties about the fitting meth¬ 
ods, especially the wavelength range. Therefore one must be 
very careful to minimize the systematic uncertainties, e.g. by 
checking the sensitivity of the results to the adopted wavelength 
ranges and other model parameters. 
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